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Abstract 

This paper discusses properties of a Doubly Stochastic Poisson Process 
(DSPP) where the intensity process belongs to a class of affine diffusions. 
For any intensity process from this class we derive an analytical expres- 
sion for probability distribution functions of the corresponding DSPP. A 
specification of our results is provided in a particular case where the in- 
tensity is given by one-dimensional Feller process and its parameters are 
estimated by Kalman filtering for high frequency transaction data. 

1 Introduction 



The D oubly Stochastic Poisson Processes (DSPPs) were introduced by ICox 



(1955) by allowing the intensity of the Poisson process to be described by a 



positive random variable and not just a deterministic function. The aim of such 
a generalization was to allow the dynamic of a process that is exogenous to 
the model to influence transitions in the point process that we are concerned 
with. Point processes have applications in several areas of applications, among 
which we mention: biostatistics, finance and reliability theory. In biostatistics 
they form a theo retical framework for studying recurrent events, as was done in 
Gail et al. (Il980f). in stu dies of the size of tumors in rats over a period of time. 



In finance. [Landol (1998) was the pioneer in using p oint processes for des c ribing 



occurrences of credit events. In reliability theory, iDalal and Mclntoshl (jl994l ) 



work has developed criteria for determining an optimum stopping time for test- 
ing and validating a software. 



It was the seminal work by [Cox] ([1955) which introduced the DSPPs (also 



known as the Cox Process) . The main work in this field is iGrandelll ([197 
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where the main properties of the DSPPs have been presented in ter ms of stan- 
dard c onstr uction of Probability Theory. Alternatively, in the books of lBremaud 
( 19721) and iDalev and Vere- Jonesl £l988) he presentation is based in line with 
the concepts and properties of Martingales. However, all these sources focus on 
deriving general properties of DSPP, without exploring the functional form of 
the intensity of the process. Consequently, they attracted a reduced number of 
applications. 



The study of the DSPPs took a new turn once the functional form for the 
intensity of the process has been specified. As a result it became possible to 
obtain analytical expressions for probability density functi ons for differen t types 



of processes. In this context we can quote may be cited iBouzas et al.1 (|2002f ) 



who made u se of truncate d norm al distribution to describe the intensity of the 
process and IBouzas et al. (2006) who generalized the form of intensity to in- 
clude the case of a harmonic oscillator. The contribution of these works was 
that closed analytical expressions for density functions of the Cox process have 
been obtained, as well as their moments. However, in both cases the authors 
used constructions in which the restriction on non-negativity for the intensity 
measure was not maintained. In order to get a round of this limitation, the 
authors defined a region in the parameter space where the probability of oc- 
currence of negative values for the intensity is reduced. With the intent to 



guara ntee the preservation of the non-negati vity condition, iBasu and Dassios 
( 2002) and Kozachenko and Pogorilva k (2008) suggested the adoption of a log- 
normal model for the process intensity. A formulati on that incorpora t es the 
intensity into a dynamic formulation is developed in iDassios and Jangl (|2008l ) 
who use the functional form of a process of the Shot-Noise type which, despite 
guaranteeing the non-negativity of intensity, is not in fact a diffusion process. 
On the other hand lWei et al.l (|2002l) assumed that the intensity is governed by a 
one-dimensional Feller process and obtained a form of probability density func- 
tion for the corresponding DSPP. 



F eller processes we re introduced and established in the financial literature 



after ICox et al.l (|1985l ). One of the properties is that the Feller process lives 
in R + which guarantees that the non-negativity condition is fulfilled. In the 
present work we assume that the intensi t y is c ontrolled by a related diffusion 
process, as formalized by Duffle and Kan ( 19961 ). which incor porates Fel l er pro - 
cesses in one or d dimen sion s. In this way, the models from | Wei et al.1 ([2002J), 
Basu and Dassiosl (|2002f ) and lKozachenko and Pogorilv ak (2008) can be seen as 



specific examples of the model proposed in the present work. 



The use of point processes in finance, especially of DSPPs, progressed con- 
siderably at the end of the 1990s with the d evelopment of models of ma nag- 
ing and pric i ng th e credit risk. In particular. iDuffie and Singleton! ( 19991 ) and 
Duffi e et al. I (120031) formalized the construction of the probability density of the 
first jump in the process, in a related context. More precisely, the time to the 
first jump represents, in the context of credit risk, the time until the bankruptcy 
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(default) of a company (and/or a country). Here, once the absorbing state had 
been reached, it was unnecessary to study the further dynamic of the DSPP. 



Recently a new area of applications of point processes in finance has emerged, 
along with the use of these models for describing the arrival process of bid and 
ask orders in an electronic trading environment. In these models, the arrival 
process of orders changes over time; the idea is to characterize the dynamic of 
this process and obtain expressions that may be treated analytically, describing 
the probability that an order had been sent in line with a market con figuration 
and was executed before the price was altered. See Cont et al~l ( 201fj) . 



Perhaps the most well known point process is the homogenous Poisson pro- 
cess. For this process the arrival rate is constant. A homogenous Poisson process 
can therefore be described by a single value Af = A. However, in many appli- 
cations the assumption of a constant arrival rate is not realistic. Indeed, in 
financial data we tend to observe bursts of trading activity followed by lulls. 
This feature becomes apparent when looking at the series of actual orders ar- 
rivals. Figure 1 presents the average number of sell orders for BRL/USD FX 
futures submitted to the Brazilian Exchange, BVMF. The plot indicates a non- 
constant behavior for the average number of submitted orders; the homogenous 
Poisson model is clearly not suitable for such data. 



Average number of sell order per minute during October 2009 
BRL/USD Futures contract at BM&FBOVESPA - Brazil 




Figure 1: Sell orders for FX futures contract during October: (ticker DOL FUT 
- expiry NOV09) 

In fact, the number of bids and asks that getting into the order book may 
depend upon a number of factors exogenous to the model. For example, it may 
be the level of investors risk aversion on a given day, or intra-day seasonality, 
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or a disclosure of a piece of information, or a new technology producing an im- 
pact on a certain sector of the economy. For this reason, in our opinion, the 
construction of a model that can be treated analytically and incorporates cn- 
dogenously a stochastic behavior of the intensity of a point process represents a 
contribution to the literature. In light of the above topics, the main aim of this 
article is to present new results on DSPPs when the intensity belongs to a family 
of so-called affine diffusions with potencial application in high frequency trading. 

In this situation seems reasonable to focus on a particular class of models and 
attempt to use the approach based on point processes in order to study their dy- 
namic over time. For this purpose we selected the Affine Term Structure (ATS) 
models, and our goal is to obtain the form Theorem for the Probability Density 
Functions for the Cox Process when the intensity belongs to a family of affine 
diffusion. Additionally, in one particular case when the intensity is governed 
by an one-dimensional Feller diffusion we obtain more detailed results, such as 
its moments and the convergence to its stationary distribution. Finally we pro- 
pose an estimation procedure for point processes with stochastic affine intensity 
based on Kalman Filter conjugated with quasi-maximum likelihood estimator. 
Thus, the estimation procedure is applied to high frequency transactional data 
from FX futures contracts traded in Brazil. 

2 The basic structure 

Consider the filtered probability space (Q,<?,G, P), where G = (Gt)t>o is a 
filtration having the sets with null measure and right continuous. 

Definition 1 Let t be a non-negative random variable t : Ci — > R + U oo 7 which 
is defined on a probability space (fl,Q,G,F), such that {r < t} e Q t Mt > 0. 

Definition 2 Given a no-decreasing sequence of stopping time {n,i <G N}, the 
Point Processes N t is a (cddldg) processes defined as: 

i 

Definition 3 Define the intensity processes as: 

d 

A* = Po + ^2pi,iX iit (2) 
i=i 

= po + pi ■ X t (3) 
With po € K and p\ € M. d . Where the state variable X t will be defined below. 

Definition 4 Let (il,Q,G,P) be a filtered probability space and let ¥ be a sub- 
filtration of G. We called cumulate intensity (or Hazard Process) an T- 

adapted, right- continuous, increasing stochastic process A ,t = J * X u du for t > 
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with P-a.s. Aq = and Aoo = oo. 



The hazard process plays an important role in the martingale approach to 
either credit risk or asset pricing because it is the compensator of the associ- 
ated doubly stochastic Poisson process, see for example Bielecki and Rutkowski 
(2002). For our purposes, the hazard process will be important to derive the 
probability density function for the Point Process, (AT t ) t > . To construct the 
stopping time with the given intensity, define as in Lando (1998): 

- : :nl {/ : A', , - /.;.,} (4) 
where is a unit exponential random variable independent of X t V i. 

Definition 5 Let X t : T> x M — > M. d , be the state variable as solution for the 
stochastic differential equation SDE: 



dX t = n(X, t)dt + a(X, t)dW t (5) 

where: (i(X, t):DxK + ^ R d , is the drift; a(X, t) : V x K+ -> M. dxd , the diffu- 
sion coefficient and, Wt the standard Brownian Motion defined in (fi, Q, G, P). 
Where T> C M. d will be properly defined below. 



Definition 6 (jDuffie and K an (1996)) Define X t as an affine processes, 

if the following condition are simultaneously satisfied for f5[): 



1. Drift 

For 9 e R d and JC G M dxd 

2. Covariance matrix 



n(X,t) =1C(G-X t ) 



(6) 



o-{X,t) = Ev^ 



(7) 



Where S is d x d matrix not necessarily symmetric, and a is a diagonal 
matrix with the z-th elements from the diagonal given by: 



1, 



,,d, 



(8) 



with aelek R d . 

Adopting the construction of Duffie and Kan ( 19961) and Dai and Singleton 



(2000) the state variable in its affine form is: 

dX t = fc(0 - X t )dt + S^dW t 



(9) 
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where Wt is the d- dimensional independent standard Brownian Motion in 
(Q, Q, G, P), JC and £ are d x d matrices which may be nondiagonal and asym- 
metric. 



According to Duffie and Kan ( 19961) the coefficient vectors 61, . . . , bd in © 



generate stochastic volatility unless they are all zero, in which case © defines 
a Gauss-Markov process. If the volatilities (all or some) are stochastic, then 
two questions face up: how to constrain the model in order to avoid negative 
volatilities? Under those constraints, which is the set Del 11 where the process 
(X t )t>o can take values? 

The open domain T> implied by nonnegative volatilities may be defined as: 



V = Lc e M. d : a, + b t ■ x > 0, ie{l,...,d}j 



(10) 



The following condition of iDuffie and Kan (|l996h is sufficient for the posi- 
tiveness of cr, : 



Condition A (Duffi e and Kanl (|l996h l: For all i: 



1. Vz such that a t {x) = 0, bj "(£(9 - X t ) > ±bjEE T bi 

2. Vj, if (bJZ)j ^ 0, then a t = kaj for k> 

Both parts of Condition A are designed to ensure strictly positive volatility, 
and they are both effectively necessary for this purpose. 

A second set of constrains must be imposed to guarantee that the non- 
negativity condition for (At)oo is fulfilled. The open domain T> implied by 
nonnegative intensities may be defined as: 

V= \xeR d :X t (x) >ol (11) 



In this sense, the strongest form of guarantee that (p~0|) and (fTT|) arc simul- 
taneously met it is to impose T> = M.~L. However, it is a very strong constraint 
on the state space resulting in a large number of models excluded. 



Using the Canonical Representation developed bv lDuffie and Sin gletonl (Il999h 



it will be possible to overcome the non-negativity for (A t ) t > . In fact, according 
to the Canonical Representation, A m (d), where m = rank(a), the state vector 
X t will be split into two subvectors: the first factors vector Xf £ R m , with 
< m < d and the other group of factors are stacked in Xf € M. d ~ m . 

This way it is also possible to overcome existence and uniqueness problems 
for a more general class of affine diffusion than the one that satisfies (|9]). In fact 
the process Xf exists and is unique, moreover it is autonomous with respect to 
Xf '. The volatility of X^ , conditionally on X t s , is given, so there is no unique- 
ness problem concerning Xf . Finally, from the Canonical Representation it is 



G 



known that the process X D , when X ^ 0, has Gaussian Distribution, inde- 
pendent from X B , with mean yP and standard deviation a D . 

From the Canonical Representation, we just need to impose an extra condi- 
tion to ensure strictly positive of (\t)t>o- 

Condition B: if we set 

fi D > , (12) 

it is possible to state that: 

Proposition 1 If^*fi D > a D , then the process X D is nonnegative almost sure, 
i.e. X D e V+. 

Proof of Proposition Q] 

According to the Canonical Representation, the process X D is Gaussian 
with mean pP and standard deviation o~ D . Thus it is possible to compute the 
probability X D eV+: 



where $(•) is the cumulative normal distribution function. 

Because condition (|12l) the problem described in (fT3|) is now written as: 



where: 7 = ^rr 

Then by requiring that 7* pP > cr D the non-negativity condition for (At)t>o 
is satisfied almost sure. 

■ 

We now turn to examples and illustrations of Afhne framework. 
• Example 1: Ornstein-Ulhenbeck (Vasicek) 




(13) 



7* = inf (7 : $(7) = 1) 



(14) 
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D 



dX t = k(6 - X t )dt + adWt 



• Example 2: Feller (or Square Root, Cox-Ingersoll-Ross) 



dX t = k(6 - X t )dt + a^XtdWt 
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• Example 3 : Geometric Brownian Motion 

dS t 



fidt + adW t 
1 



St 



X t := lnS t dX t = Ut - |cr 2 ) dt + adW t 

Example 4: Heston (Stochastic Volatility) 



X t := lnS t dX t = [n - -a 2 ) dt + a t dW} 



do t — k(6 — o t )dt + v^faldW 2 

< dW t \dW? >= p 
• Example 5 : Multivariate CIR 



dX t = 





K\6l 




' -Ki 


( 


K2&2 


+ 


-K 2 



X t dt- 



ai 
a 2 



1 
1 



X t dW t 



In order to have the no-negative random variable 7$ as a stopping time 
a technical conditions must be imposed over the filtration. So we shall first 
examine a trivial case where the condition is not respected. 

Proposition 2 The hitting time t is not a F-stopping time with respect to 
T t := a{W s : < s < t} 

Proof of proposition [2] 

Assume, by absurd, that r is F-stopping time for Tt '■— cr{W s : < s < f}, 
so using the Martingale Representation Theorem there exist a compensated 
process (a F-Martingale) M t which may be represented as a stochastic integral 
with respect to the Brownian motion W . Therefore we arrive a contradiction 
because M t must jump in r. 



Thus the filtration F should be enlarge. There are several ways to expand it 
but we do enlarged F just to get r as stopping time. Then the proper filtration 
G = {Gt)t>Q will be constructed as: 



wher43 



Gt=TtV H t 
Ut = o{N s :0<s<t} 



A\J B:= maximum between A and B 
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Definition 7 Let (Q,G,G, P) be a filtered probability space and N t a stochastic 
processes defined on it. Let ¥ be a sub-filtration ofG, then N t is called a Doubly 
Stochastic Poisson Process with relation to A t , if At is ¥ -measurable and for all 
< t < T and k = 0, 1, . . . the next condition be satisfied: 



P{N T -N t = k\U t V7 co ) = ^ (A t , T ) fc exp (-A t , T ) 



(15) 



where Foo — a(T u : u G M + ) 
As a particular case: 

P{N T -N t = k\U t V Too) = P(N T -N t = k^) 



(16) 



Therefore, conditioned to cr-field Foo the increments of Nt — Nt are inde- 
pendent of the cr-field Ht- 

Thus, the probability of no occurrences within the interval [t,T] for the 
processes N t with intensity A t > is given by: 



P{N T - N t = 0) = Et 



exp 



X„du 



(17) 



2.1 Laplace transform for the process A t 



The Laplace transform of stochastic processes is a key ingredient to achieve our 
results. However Laplace transform for integral of stochastic processes may be 
obtained in closed form only for a limited number of processes. 



Alb anese and Lawil (120041 ) formalize the criteria to define which processes 



have a analytic form for its Laplace transform. Thus take a diffusion (Xt)t>o 
defined on (fl,Q,G, P) and consider the Laplace transform defined by: 



L t (ji,X,t) :=Et 
where t < T, fi g C and / 



f{X, s)exp y-Vj <j){X,s)ds ^ 
(4:MkM two Borel functions. 



(18) 



Thus, it is possible to state that: 

Result 1 ( Alba nese and Lawil ([2004)) The class of stochastic processes with 
Laplace transform for its integral is given by: 



dX, 



t h'(X t )A(X t ) 2 dt i V2A(X t ) 



h(X t ) R{X t ) 



VWXt) 



dW t 



(19) 



With additional conditions: 
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1. Three second order poly onymous in x: A{x),R{x),Q{x). Such that A(x) 
belongs to the set o {l,x,x(l — x),x 2 + 1} e R(X t ) > 0; 

2. The function h(x) is a linear combination of hypergeometric functions of 
the confluent^ 

From type \F\ if A(x) S {l,x} and gaussian from type 2-F1 otherwise. 
Thus the Laplace Transform L t (fi, X,t) is defined by: 

^ = jR(xJ (20) 



Example 6 : As an application from the above result we shall prove that the 
Laplace transform of A f :— J Q X u du when the intensity follows a one-dimension 
Feller process exists. 



Thus assume that the polynomials are defined as: 



A(x) = x, R(x) = 



2x 



h(x) = x 



(21) 



Substituting the polynomial into (|19|) with a variable change we have: 
d\ t = k{6 - X t )dt + o-^/\~tdWt 

So in those cases where the Laplace transform L t (/J,, X, t) exists it is possible 
to apply the well know resullH: 

Result 2 (Feynman-Kac) Let (X t ) t >n be a diffusion process with infinitesi- 
mal generator A. Assume that F £ C 2,1 (T> d x [0,t)) and V £ C 1 is bounded. 
Then: 



f(Xt) exp — fi / V(x, s)ds 



F(x,t,/J,) 



it is solution of the partial differential equation (PDF): 



dF 
dt 

F(x,0,n) 



AF(x, t, fj.) - V(x, t)F(x, t, n)) 

fix) xev 



(22) 



(23) 



2 A hypergeometric function in its general form may be written as: 
p F q {a 1 ,. . .,a p ;ji, ...,-y q ;z) 
For p < q -f- 1, 7j £ C \ Z+ been represented using Taylor's expansion around z = 

(c*i)„ . . . (q p )„ z n 



pF q (a!, . . . ,a p ;7i, 



.7« 



;*) = £- 



V (Ti)n • • ■ {y q )n n! 
3 Details might be found in Karatzas and Shrcvc ( 1991) 
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and 



dF dF 1 
AF( X ,t) = - + -,( X ,t) + -Tr 



a(x, t)a(x, t) 



d 2 F(t,x) 



(24) 



dXdX 

Furthermore 0, ifU{x,t,n) £ C 2 ' 1 solve fS3j), then U(x,t, p) = F(x,t, p). 

Taking the Theorem [2] it is possible to state that: 

Proposition 3 Let (X t )t>o be a diffusion processes satisfying the regularity 
condition so the solution to (23\) has the form: 



F(x,t) = e °(*)+/s(t)-* (25) 
where the coefficients a(t) and (3(t) are deterministic and satisfying: 



fi'(t)=p 1 -K T m-^ T ^ (26) 
a\t)=p Q -1Cep{t)-\p T aP (27) 

with a £ M and b £ R d . 
Proof of Proposition [3] 

Imposing that the processes (X t )t>o and (\t)t>o are afhne as © and ([2]) so: 

o = -(p + Pv x)F(x,t)+F t (x,t)+F x (x,t)(ic(e-x t ))+± g d Qxf d x- fa+kj-x) 

(28) 

Inserting F(x,t) — gP'W+Pw-^ m to the PDE above and grouping the terms 
in x: 

u(-)x + v(-) = 

Where 

u(.) - -/3 1 (t) +px- JC T (3(t) - ^/3(i) T 6/3(t) (29) 

= a 1 (t) + a, - £90(4) - i/?(i) T a/3 (t) (30) 

Use the separation of variable technique to obtain that a and /3 satisfy a Riccati 
equation with boundary condition a(0) = 0e /3(0) = w. 



4 We write Tr for trace of a matrix. 

5 where 1 stand for the derivative with respect to t. 
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3 Distribution of Doubly Stochastic Poisson Pro- 
cess 

For a non-homogeneous Poisson processes (N t )t>o with intensity (\t)t>o, the 
probability of k occurrences within the interval [t, T] is given by: 



P(N T -N t = k) = ~E 
k\ 



T y / fT 

X u du cxp — / \ u du 



(31) 



for k = 0,1,2, ... 

Thus it is possible to state the important result: 

Theorem 1 The Probability Distribution Function for a Doubly Stochastic Pois- 
son process, (N t )t>o> within at interval [t,T] when the intensity is an affine 
diffusion as (0|) is given by: 



P(N T — N t = k) = — E 
k\ 



(A t ) fc e 



fc! 



(32) 



where G\ (//) is the k-th derivative of Moment Generating Function for the 
Hazard Process, At 

Proof of Theorem [1] 

Imposing f(X t ) = 1 and V(X t ) = A t on the left hand side of equation (|22[) 
it may be seen as the Laplace Transform (or Moment Generating Function) for 
the Hazard Process, At := J X u du. 

Thus, using that if G x (p) is the MGF of X, then 



d k G x {n) 



d/i k 



G k M=^{X k e» x ) , 



(33) 



and so the result follows. 



In a particular case when the intensity follows an one-dimensional Feller 
process we have: 

Theorem 2 The Probability Density Function (PDF) for a non-homogeneous 
Poisson processes (N t )t>o within the interval [t,T] when the intensity takes the 
form 

d\ t = k(6 - \ t )dt + o-yf\~tdWt (34) 

is expressed by: 

P(N T - N t = k) = {A t ) k e 



k\ 



fc! 



G At (l) 



(35) 
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G At (1) = E (e" A T A "* 1 ) = e «(*.T)-^(t,T)A t 



alt, T) = — 5- m 



27 (e 



(7+«0/2^| 



(7 + k) (e-r( T -*) -1) +2 7/ 
2/^(e-^ T - t )-l) 



(7 + k) (e-^- 1 ) - 1) + 2 7 



where 7 = \/ k 2 + 2a 2 



Proof of Proposition [2] 

See appendix A 



(36) 

(37) 

(38) 
(39) 



Since DSPPs processes are essentially Poisson processes, each result for Pois- 
son processes generally has a counterpart for DSPPs processes. The following 
are some basic pro perties of a DSPP process ( Nt) f .>n driven by (A t ) t >o- See 



are some basic pro perties or a usff process ( i\t)t>n driven oy (A t )t 
iBremaudl (fl972h or iDalev and Vere- Jonesl (|l988h for further properties 



Proposition 4 The first two moments of (N t )t>o where the intensity is an 
one- dimensional Feller process are: 



1. E(JV t ) 



et 



(40) 



2.Var{N t ) = 

(e- Kt +l)(\ -8)-2(8+\ ) 
Proof of Proposition [4] 



26t r 

K 



a 2 fOe 



4e" K - 5 



A e 



-2k* 



The first moment is obtained using the Laplace transform of A t : 
For the second moment we need an additional result: 



Var (A* 



d_ 

fi=0 \ BfJL 



oH /30-2A o 



(41) 



-f/o A = ds 



/x=0 



(42) 



Substituting in Var(iV t ) = E(A t ) + Var(A t ) we obtain the result. 
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The stochastic nature of the intensity causes the variance of the process to 
be greater than the variance of a homogeneous Poisson process with the same 
expected intensity measure. This feature of the DSPPs processes is referred, in 
the literature on point processes and survival models, as overdispersion. 

4 Stationary distribution for (Nt)t>a 

Stationarity is a very important concept in time series analysis. A stationarity 
assumption will allow us to estimate parameters from point processes and make 
predictions. The characterization of stationarity for a Point process relying on 
whether its intensity process is stationary. 

Definition 8 A point process (N t )t>o is stationary (or isotropic) if for all 
Ai, . . .,A„ eJ e/iel 

(N Al+h ,...,N An+h )±(N Al ,...,N An ) 

Thus, the increments of N t are translation invariant distribution with respect 
to any translation h G R. 

The definition above can be written in a short form: 

9 h N = n, hem. 

Where: 

e h N A ■= N A+h 

and = means equally in distribution. 



Proposition 5 The Doubly Stochastic Poisson Process {N t )t>o is stationary if 
its intensity (Xt)t>o * s stationary. 

Proof of Proposition [5] 

The Laplace transform for a point process, (N t ) t >o, with intensity process 
(A t ) t >o is given by: 



L N (f) =E exp 



(l-e-W>)\(t)(dt) 



(43) 
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Thus, taking a/iel with (|35]l: 



igfcjv = IE |exp 
= E jexp 
= E < cxp 



- / f(t)N(dt + h) 



- I f(t-h)N(dt) 



- / (1-e 



-f(t-h) 



)X(t)(dt) 



= e|cx P - J^l-e-f^)6 h \{t){dt) | (44) 



Therefore, if (Xt)t>o is stationary such that 6hX — A, so together with (|44l) we 
obtain: 

L ehN (f) = L N (f) (45) 



In a recent paper IGlasserman and Kiml (|2010| ) analyze the tail behavior, the 
range of finite exponential moments, and the convergence t o stationarity in affine 
model s, focusing on th e class of canonical models def ined bv lDuffie and Singleton! 
( 1999f h According to IGlasserman and Kim ( l20ld) ,he one-dimensional Feller 
process has a stationary distribution so the next step is determinate it. 

Proposition 6 The stationary distribution of one- dimensional Feller process 
is the Gamma distribution with parameters 

Proof of proposition [6] 

See Appendix B 



Definition 9 Suppose that u, fj,%, . . . are locally finite measures defined on (f2, Q, G, P), 
a necessary condition for Vague Convergence of [i n to a, in short /i n A a, 
is 

IJ-nf -> fJ-f, when n -> oo, for f £ C£(fi) 
Where: C£(fl) is a continuous function with compact support / : O — > K+ 

From the definition |H] we shall present without proof the well know result for 
convergence of point process: 

Result 3 Let N\, N2, ■ ■ ■ defined on (f2, Q, G, P) be a sequence of Point processes 
the following results are equivalents: 

1. Nt A 



2. N t f A N^f, for all f e 
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3. E (e N 'f) A E (e™^) , for all f £ C+(fi) 
Proof of Theor em [3] 



Kallenberd |l986) 



Theorem 3 For all t > 1, suppose that (Nt)t>o is Poisson process defined 
on (SI, 0,G,P) twi/j intensity (Xt)t>o- If At A A and A is locally finite, thus 
N t A A^oo, where Noo is a Poisson process with intensity A 

Proof of Theorem H 

The Laplace transform for a point processes (02]) may be written as E [e iVt -'] = 
e Xth , for / e C£(0), where fc(t) = (1 - e _/ W). The Theorem 1 establish that 
A t /i A Aft-, so: 

E [e^'] = A e Ah = E [VW] (46) 
It follows from result [3] that Nt A JVqo . 



From the above results it is possible to determine in a closed form the sta- 
tionary distribution of (N t )t>o- We now turn to one application of this result. 

Proposition 7 The stationary distribution of (N t )t>o with intensity (Xt)t>o 
given by dXt = k(6 — Xt)dt + o^/XtdWt is the negative Binomial distribution. 

Proof of Proposition [7] 

The proof is straightforward when Proposition [5] is considered together with 
Result 1 



(N t =k)= F(N t = k\A = X)^{X)dX 







ri^^-e-^dX 

o k\ T(a) 

t k F* T(a + k) f°° (!3 + tr+ k x(0+t) k+a _ 
k\T(a)(l3 + t) a + k J T(a + k) 



*(N t = k) 



(a + k- 1)! / t \ k ( (3 



(a-l)!ar! \P + tJ \P + t 
Afeg(a,p) with p = -—-^ 



16 



We have determined the stationary distribution for (N t )t>o when the inten- 
sity is given by an one-dimensional Feller process. Closely related to this result 
is the determination of the rate of convergence of (N t )t>o to its stationary dis- 
tribution N£ . 

Theorem 4 The Cox Process (N t )t>o converge to its stationary distribution, 
N*, exponentially fast at rate 2k. 

Proof of Theorem [4] It is known that the transition density for the process 
(^t)t>o when it follows a Feller process is a non-central ^-distribution. From 
proposition [6] we have sho wn that the proces s (\t)t >n converge to the Gamma 



distribution. According to iKarlin and Tavlorl (|1981f ). it is possible rewrite both 



distribution by its spectral representation: 



p{t,x,y) - (2 K ) y e L n [ a2/2 )K { a 2 /2 ) T (n + k9) 



(47) 

where L„ K is the Laguerre polynomial with parameter 2k — 1, and with 
the following property Lq K_1 (-) = 1. T is the Gamma Function. 

To describe the stationary distribution p(y) by its spectral representation 
we just need to set n = at equation (|4T)l. 

We decided to evaluate the convergence speed of (A t ) t >o by looking at the 
convergence of its density functions written in the spectral form: 

\p(t,x,y) -p(y)\ = <D(e- 2Kt ),\Jx > when t -> oo (48) 

Because according to the spectral representation, we have that for Vn, the 
only term inside the summation involving t is e~ 2Kt . 

Finally, the unconditional version of Cox process gives us: 



— (A t)T ) 'exp(-A t!T )p(t,a;,?/)|J : ' t 



(49) 



P(N T -N t = k) = J 
Therefore, 

\N t -N?\ = 0(e- 2Kt ) when t -> oo (50) 



5 Model Application 

In section 2 we developed the requisite theory used to construct the Cox Process 
as an afhne function of the underlying state variables. In every case, this rela- 
tionship was subject to a given parameter set. Unfortunately, the theory does 
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not tell us anything about the appropriate values that must be specified for this 
parameter set. We must, therefore, turn to the econometric literature to handle 
this important issue. Since the seminal paper by Engle and Russell (1998) the 
modelling of financial point process is an ongoing topic in the area of financial 
econometrics. The financial point processes are associated with the random 
arrival of specific financial trading events, such as transactions, quote updates, 
limit orders or price changes observable based on financial high-frequency data. 

Moreover, it has been realized that the timing of trading events, such as the 
arrival of particular orders and trades, and the frequency in which the latter 
occur have information value for the state of the market and play an important 
role in market microstructure analysis. 

Although the literature on the parametric estimation of point processes (fi- 
nancial Point Process as well) is as large as the theoretical literature, there is 
as yet no consensus as to the best approach. Based on that, we propose a new 
technique to the estimation of Cox Processes parameters. The methodology we 
will be using, based on the Kalman filter, exploits the theoretical affine relation- 
ship between the Cox Process and the state variables to subsequently estimate 
the parameter set. The strength of this approach is that Kalman filter is an 
algorithm that acts to identify the underlying, and unobserved, state variables 
that govern the Cox Process dynamics. Once the unobserved component is fil- 
tered the Quasi-Maximum likelihood estimator will be able to determine the 
model parameters. 

In order to estimate parameters and to extract the unobservable state vari- 
ables we restrict the equation ([35]) to deal with the probability of no arrivals, 
P(Nt — Nt = 0), within the interval (t — T). Therefore, the probability of no 
arrivals for the Cox Process with Feller diffusion is 

P(N T -N t = 0) = e «(*,T)-/J(t,T)A. (5i) 
where a(t, T) and /3(t,T) were defined, respectively, in equations ([57)1 and (|55|) . 

Additionally it is possible linearize (|51l) : 

In P(N t , T = 0) = In a(t, T) - f3(t, T)X S (52) 
Thus, the measure equation is log-linear in At and it can be written as: 

\nP(N t , T = 0) = \na(t,T) - /3(t,T)\ s + Xs (53) 

where: 

Xs~N(0,TL s ) 

The inclusion of an error term in equation (|53[) is motivated by the fact that 
the underlying intensity process may be inadequate. If the true factor process 
is not a Feller process equation (|5Tj) will be functionally misspecified and esti- 
mates of At will be inferior. In this case the probability of no arrivals within 
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the interval (t, T) implied by the Feller process will systematically deviate from 
observed arrivals. Therefore, in a correctly specified model the errors Xs should 
be serially and cross-sectionally uncorrelated with mean zero. 

It is known that the exact transition density p(t,x,y) — P(X t € dy\X t -i = 
dx) for the Feller Process is the product of K non-central x 2 -densities. Esti- 
mation of the unobservable state variables with an approximate Kalman filter 
in combination with quasi-maximum-likclihood (QML) estimation of the model 
parameters can be carried out by substituting the exact transition density by a 
normal density: 

A S |A S _ X ~ N(p„Q.) 

where fi s and Q s are defined in such a way that the first two moments of the 
approximate normal and the exact transition density are equal. The moments 
are time varying and defined as: 



Us = 9[1 — exp(— k)] + exp(— k)A s _i 
and Q s is diagonal matrix with elements: 

2 1 - exp(-K) 



-[1 — exp(— k)] + exp(— k)A s _i 



(54) 



(55) 



Note that once we have represented the C ox Proc e ss wi th Feller intensity 
into a State Space form we can, according to lHarvev to calculate the 

quasi log-likelihood: 



T T 

log/(y t |x t ;0) = -ilog2 7 r(T-if)-i ]T log|S t |-^ £ y' t S^y t (56) 

t—K+l t=K+l 

where T is the sample size; K is the dimensional of the Feller process. 



6 Empirical Analysis 
6.1 Data description 

In this section, we apply our estimation technique to some transaction data 
from the Brazilian Exchange (BM&FBOVESPAjj. 

The sample is formed by all submitted sell orders for BRL/USD FX futures 
contraciQ traded in BM&FBOVESPA during October 2008 with expiry date of 

6 BM&FBOVESPA is the fourth largest exchange in the word in terms of market capital- 
ization. BM&FBOVESPA has a vertically integrated business model with a trade platform 
and clearing for equities, derivatives and cash market for currency, government and private 
bonds. 

7 Ticker: FUT DOLX08 



19 



November 1, 2008. This FX contract is one of the most liquid FX contracts in 
the emerging markets and the average volume of 300.000 traded daily is signif- 
icant even for developed markets. We have a total of 535 records with trading 
occurring continuously from 10 am to 6 pm exclusively through the electronic 
venue. 



The BM&FBOVESPA electronic trade system (GTS) uses the concept of 
limit order book, matching orders by price/time priority. Lower offer price take 
precedence over higher offers prices, and higher bid price take precedence over 
lower bid prices. If there is more than one bid or offer at same price, earlier bids 
and offers take precedence over later bids and offers. The offers are recorded 
in milliseconds allowing the highest precision for determination of precedence 
criteria. We also observe that no two consecutive orders arrive in the order 
book in an intervals smaller than 10 milliseconds, probably due to the internal 
network latency 



While the probability of no arrivals are not themselves directly observable, 
we can use the empirical frequency of orders arrivals as proxy. To construct the 
empirical frequency we aggregate orders sent within a same minute for every day 
during October and the average value is taken. Thus, due to the 10 milliseconds 
network latency, it is possible, at least theoretically, that 6.000 sells orders could 
arrive in the order book in 60 seconds. 

Therefore, the observable variable is constructed at@: 



2/(t- s =60 sec) 



#(Scll orders) 
6.000 



(57) 



6.2 Estimation results 

In section we illustrate the properties of the Kalman filter for estimating the 
parameters of Cox Process with Feller intensity. Table 1 contains the parameter 
estimates = (8, k, a) . The estimated standard deviations of errors - the square 
root of the diagonal elements of H - are also presented in Table 1. Standard 
errors for the QML estimates are obtained as described in Hamilton (1994, 
p.389). 





6 


K 


a 


Std. Error x 


Estimates 
Std Error 


0.065 
2.6E-05 


0.0043 
1.7E-05 


0.00267 
7.47E-08 


0.0010 
2.68E-09 



Table 1: Estimated parameters of the Cox Process with affine intensity applied 
for BRL/USD Futures sell orders data: (ticker DOL FUT - maturity NOV09) 



where # is the count measure within the interval (t — s = 60sec) 



20 



We obtain highly significant parameter estimates (at the 1% level). Addi- 
tionally, the estimated standard deviation for the measurement error x is twice 
smaller than the diffusion parameter a. 



x - estimated o- actual 
2.2 1 1 1 1 r- 




100 200 300 400 500 SOO 



Figure 2: Cox Process with Feller intensity fitted to BRL/USD futures sell 
orders (ticker DOL FUT - maturity NOV09) 

From figure 1 we can see the model flexibility in reacting to different changes 
in the numbers of orders submitted. In this case, the model seems to fit quite 
well to actual order flow and the differences between actual and theoretical 
probability indicate that, on average, the model tends to slightly under-estimate 
actual occurrence of sell orders. In order to assess the model fitted we conducted 
a diagnostic checking for possible misspecifications based on the standardized 
residuals. For a well specified model we need the residuals follows a white noise 
process. 
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L 


Q(L) 


p- value 


5 


8.36 


0.13 


10 


13.47 


0.19 


15 


17 


0.31 



Table 2: Ljung-Box test for residual correlation 



Table 2 contains the Ljung-Bo?H statistics up to the 15" 1 lag. The LB test 
statistic allows us to not reject, at a high significance level, the hypothesis that 
errors are white noise. In the absence of dependence in errors the Feller process 
for intensity seems to be a right choice to describe the observed sell orders 
behavior. 

7 Simulation results 

In this section, to analyze the performance of our estimation algorithm we simu- 
late various Cox Process outcomes using a known parameter set and proceed to 
estimate the model parameters. This simulation exercise is intended to indicate 
how effective this technique is in terms of identifying parameters. The simula- 
tions for the unobservable state variables have been drawn from the noncentral 
Qui-Square distributiorFl. f\(d, I), and the measurement errors have been sim- 
ulated as normal random variables with zero means. 

The simulation of a sample path for the Cox process with 500 elements, 
followed by application of the estimation algorithm, is repeated 250 times. The 
following table summarize the results of the simulation exercise for the Cox 
Process. The table reports the true values, the mean estimate over the 250 
simulations, the Mean Quadratic Error, and the associated standard deviation 
of the estimates. 

9 Thc Ljung-Box statistic tests the hypothesis that a process is serially uncorrelated. Under 
the null hypothesis, the Q(L) statistic follows a Qui-Square distribution with L degrees of 
freedom, where L is the maximum number of temporal lags. 

10 The samples were drawn from the Chi-Square distribution with d degrees of freedom and 
noncentrality parameter, I: 

_ a 2 (l -e-<*s)) 



I 



(7- 

Ase _ K ( t _ s) 
Cs 
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9 


K 


a 


Actual Value 


0.04 


0.2 


0.05 


Mean estimate 


0.04 


0.3698 


0.0448 


MQE 


1.5303e-005 


0.37 


0.0023 


Std. Error 


0.0039 


0.52 


0.0001 



Table 3: A Analysis of the Kalman Filter estimates by Monte Carlo Simulation 



Distribuigao de 8 
0.8 1 1 i 1 1 




Figure 3: Empirical distribution of 9 



Distribuigao de o 




a I — < 1 ' 1 1 1 

0.04 0.041 042 043 044 045 046 047 048 

Figure 4: Empirical distribution of a 
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Distribuigao de k 

0.015 I 1 1 1 1 — 



01 



CO: 








Figure 5: Empirical distribution of k 
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8 Final Remarks 



Making use of the theoretical framework developed to model interest rate term 
structure allow us to obtain the Probability Density Function for a Double 
Stochastic Poisson Process (DSPP) when the intensity process belong to a fam- 
ily of affinc diffusion. Furthermore, the stationary distribution for DSSP may be 
found whenever the intensity process is also stationary. To illustrate our results 
in this paper one of most common diffusion in interest rate dynamics literature 
is assumed to drive the intensity process, the one-dimensional Feller process. 
However the results derived here are valid for any type of affine diffusion in 
d-dimension. 



Finally this paper does not have an empirical focus, and these results are 
primarily illustrative. The empirical analysis was included to highlight that 
applied papers on estimating the parameter set for the models examined in 
section 2 are straightforward. 



A Proof of proposition [16 

Specifying the form of /i(x,t) and a(x,t) into the Theorem [2] we have: 
dF ,„ „dF 1 7 d 2 F , 

With boundary condition F(T, A, /i) = 1 



(58) 



Using the results from proposition [3] together with (|58|l"1: 



—BX 



F x = -ABe 
F\\ = AB 2 e- BX 
F t = A t e- BX - AB t \e- BX 



Plugging the above results into 

r 2 



A (y Ab2 ~ ABt = k6AB ~ A t ( 59 ) 

Since the left hand side is a function of A, while the right is independent of it, 
the following equations must be satisfied: 

^-B 2 - B t + B - 1 = (6Q) 
A t - k9AB = { ' 



11 Where subscript denotes partial derivative 
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The first term of (|60p is Riccati equation with solution B(t,T) = v(t, T) /u(t, T) 



where v(t,T) e u(t,T) are solutions to the following system: 

f ^v(t,T) + u'(t,T) 

I u{t,T) + v(t,T) - Kv(t,T) = 



Let A = T — t, so -St = — -Ar and the system above may be written as: 



I ^v(A) - u(A) 

1 u(A) - v'(A) - kv(A) 

From the second term of (l62l) we have: 



u(A) = + v (A) + Ku(A) 
u'(A) = + u"(A) + «r/(A) 

Substituting into (|62]l and rewriting this in terms of P-operators: 

~2 ~ 



£H + kD ■ 



v(A) = 



(61) 



(62) 



(63) 



(64) 



Taking the roots of the quadratic equation (f65|) together with the boundary 
condition allow write the solution as: 

v(A) = e °' 5 (T-«) A - e o.5(- 7 -«)A (65) 

Substituting into (j6"3")) gives: 

u(A) = 0.5(7 - K)e a5(7 ~ K)A - 0.5(- 7 + K )e°- 5 ^- K > A (66) 

Since A = T — t, the solution of the Riccati equation is obtained from (|65p and 
(55): 



B(t,T) = v(A)/u(A) 

2 ( e -7(T-t) _ 

B(t, T) - 



( 7 + k) (e-^( T "*) - 1) + 2 7 



(67) 



Now consider equation (JBDJ) with T fixed, so A(t, T) is a function of i only. 
Hence: 

^ = k9AB 
dt 



A(t, T) = exp (^- K ° J t B ( s > T)ds^ 



inserting B(t,T) according to (|67)l gives: 

12 ?i e v are functions of t and T, but T is fixed, hence v (t, T) denote the derivative with 
respect to t 
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A(t, T) 



2k6 



In 



2 7 ( e (7+«)/2) 

( 7 + k) (e^( T -*) - 1) + 2 7 



(68) 



B Sketch of the proof of proposition [6] 

Let p(s, t, x, y) be the transition density of {At : t € for simplicity assume 

that X t = c E R + at instant 0. From the Forward Kolmogorov equation: 



dp 



i d 



2 r 



(69) 



The stationary distribution must satisfies ^(y) — J ^(x)p(t, x,y)dx for all t> 

dp. 



0. When reached the stationary distribution is independent of i, so 757 (i, y) = 0. 



Combining it together with (|6"9")l . we havj 13 !: 
1 d 2 







2<9y 2 



a\y)^{y) 



(70) 



Let s(y) = exp ^— dvj be the integrating factor, so integrating ([70]) 



twice: 



*(x) = C; 



c 2 



f 



s(x)a 2 (x) s(x)a 2 (x) 
m(x)[C 1 S(x) + C 2 ] 



(71) 
(72) 



Where: S(x) = J x s(y)dy 



The constants C\ e C2 are determined such that ^>(x) be a probability 
density function, so 



1 



1 



Mcr 2 (cc)s(cc) 



c?x V 16 (Z, r) 



Thus, using the form of s(x) and simplifying we have: 

1 



¥(A) = /3 C 



T(a) 



In short, ^(A) = Gamma(a, 



3 A rigorous proof might be found in iPinskvl jl995t l p. 181 and 219 
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